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UNSTEADY EULER EQUATIONS 


In transonic flutter problems where shock motion plays an important part, 
it is believed that accurate predictions of the flutter boundaries will require 
the use of codes based on the Euler equations. Only Euler codes can obtain the 
correct shock location and shock strength, and the crucially important shock 
excursion amplitude and^ phase lag. (For a discussion of the importance of 
shocks in transonic flutter, see Ref. 1.) The present study is based on the 
finite volume scheme developed by Jameson and Venkatakrishnan (Refs. 2,3) for 
the two-dimensional unsteady Euler equations. The equations are solved in 
integral form on a moving mesh, Eqs . (1-2). Here the variables p, p, u, v and 
e are the pressure, density, cartesian velocity components, and total energy, 
respectively, and and y^ are the velocity components of the moving boun- 
dary 3ft of an element ft , By applying Eq. (1) to each element or cell 
(i,j), a system of ordinary differential equations is obtained, Eqs. (3), where 
Sjj is the cell area, Q^j is the net flux out of the cell, and Dy repre- 
sents dissipative terms added to damp numerical oscillations (see Refs. 3,4). A 
five-stage Runge-Kutta scheme is used to integrate Eqs. (3) forward in time. 


FINITE VOLUME FORMULATION 


• INTEGRAL FORM ON A MOVING MESH 
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TYPICAL SECTION MODEL 


The wing is modeled as a typical section, with two degrees of freedom 
(bending h and torsion a), as illustrated in Fig. 1. The usefulness of this 
model in capturing the fundamental features of bending-torsion flutter is by now 
well established. In the usual notation, the equations of motion are of the 
form given by Eqs. (5) and (6), where the lift and moment coefficients Cl and 
depend on the motion of the airfoil. Because we will consider finite 
(rather than infinitesimal) amplitude motion, the superposition principle cannot 
be used. In the present study, Cl and C^ are calculated numerically from 
the unsteady pressure coefficient on the airfoil surface at the end of each time 
interval, obtained from the numerical solution of the Euler equations. It 
should be emphasized that the equations of motion are nonlinear through the 
dependence of Cl and C^ on the motion h,a of the airfoil (and its time 
history) . 
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METHOD OF SOLUTION 

Aeroelastic stability is determined by integrating the equations of motion 
for the coupled fluid-structure system. The structural equations are first 
transformed to normal coordinates r) r • Eqs • (7)-(8), where the columns of [$] 
are the eigenvectors of the free vibration problem. The structural integrator 
is based on the convolution integral solution, Eq. (9), and the generalized 
aerodynamic forces Q r are assumed to vary linearly within each time step At . 
Because the multi-stage Runge-Kutta scheme used to integrate the unsteady Euler 
equations was found to be sensitive to the manner in which the airfoil boundary 
condition was updated and the mesh moved, the structural integrator has been 
imbedded within the Runge-Kutta scheme in the Euler code. This permits an ef- 
ficient implementation of the exact airfoil boundary condition, Eq. (10), on 
the instantaneous position of the airfoil, given by B(x,y,t) = 0. Nonref lective 
boundary conditions are used in the far field. 


• COUPLED EQUATIONS FOR FLUID & STRUCTURE ARE 

INTEGRATED NUMERICALLY USING NORMAL COORDINATES 
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STRUCTURAL INTEGRATOR IS BASED ON CONVOLUTION INTEGRAL SOLUTION 
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• STRUCTURAL INTEGRATOR IS IMBEDDED IN FIVE-STAGE 
RUNGE-KUTTA SCHEME FOR EULER EQUATIONS 

- EXACT AIRFOIL B.C. IS SATISFIED 
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- MESH IS MOVED AT EACH STEP 
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NUMERICAL RESULTS 


Flutter calculations have been carried out for the three aeroelastic test 
cases listed in Table 1, and compared to previously published calculations based 
on various TSD codes. Case A is the same as studied by Isogai (Refs. 5,6) and 
later by Edwards et al. (Ref. 7) and also by Weatherill and Ehlers (Ref. 8). 
Note that the elastic axis location "a" is ahead of the leading edge; the idea 
here is to simulate the vibratory behavior (in pitch and plunge) of the stream- 
wise sections near the tip of a swept-back wing. Case B has been studied pre- 
viously by Isogai (Ref. 6) and by Ueda and Dowell (Refs. 9,10). Case C was 
introduced by Ueda and Dowell as an example where nonlinear (amplitude) effects 
were clearly discernible, based on LTRAN2 aerodynamics implemented via the 
describing function method. In all cases, the airfoil is forced for 3-6 cycles 
in pure torsion at a reduced frequency of interest, released, and the aero- 
elastic equations are integrated forward in time for another 3-6 cycles. The 
flutter boundary is located by calculating the logarithmic decrement 6 of the 
transient solutions, and interpolating to 6=0 between adjacent solutions with 
different U/b« a . 


TABLE 1 


Aeroelastic Test Cases 


Case 

A 

B 

c 

Airfoil (s ) 

NACA 64A010 

NACA 64A010 

NACA 64A006 

M 

0 

1 

o 

0.80 

0.86 

00 

a 

-2.0 

-0.3 

CO 

o 

1 

x a 

1.8 

0.5 

0.5 

r 2 

3.48 

0.49 

0.49 

a 


60 

60 

60 

V"a 

1.0 

0.2 

0.2 

Refs . 

5, 6, 7, 8 

6,9 

10 
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MESH GENERATION 


The unsteady Euler calculations are carried out on a C -mesh of quadrila- 
teral elements, generated by means of a square root transformation followed by 
selective stretching to compress the grid near the trailing edge. A near field 
view of the resulting mesh is shown in Fig. 2. In the far field, the mesh 
extends to 15-100 chords, depending on direction. The mesh moves with the air- 
foil as a rigid body, i.e. without deformation. Flutter calculations published 
earlier by the authors (Ref. 4) were carried out on a 96 x 16 C-mesh, which 
was found to give adequate engineering accuracy in most, but not all, of the 
cases studied. In the present study, additional calculations have been performed 
on both 96 x 16 and 192 x 32 C-meshes, and the results of Ref. 4 have been 
updated where appropriate. 



FIGURE 2 
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FLUTTER BOUNDARIES FOR CASE A 


Previous studies (Refs. 58), which have been based on a number of dif 
ferent transonic small disturbance (TSD) theories, generally agree that the 
flutter boundary for Case A exhibits a significant "transonic dip", as shown in 
Fig. 3. Also shown in this figure are the results of flutter calculations based 
on the present Euler code, and using an initial forcing amplitude of 0.1 degree 
in pitch. Overall, the agreement with previous TSD calculations are fairly 
good. However, the Euler calculations appear to shift, the bottom of the 
"bucket" toward higher Mach numbers. It is interesting to note that the bend- 
back of the flutter boundary around M - 0.88 observed by Edwards et al . (Ref. 
6) and Weatherill and Ehlers (Ref. 8), is also predicted by the present Euler 
calculations. Not surprisingly, the precise location of the nose of the curve, 
where the flutter boundary has a vertical tangent, was found to be sensitive to 
the mesh size used in the calculations. 



MACH NUMBER 


• Present code; 96 x 16 mesh 

# Present code; 192 x 32 mesh 
O HYTRAN2 (Ref. 7) 

V 0PTRAN2 (Ref. 8) 

A EXTRAN2 (Ref. 6) 


FIGURE 3 
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FLUTTER FREQUENCIES FOR CASE A 


Figure 4 shows the flutter frequencies vs. Mach number, with some com- 
parisons to earlier 0PTRAN2 calculations by Weatherill and Ehlers (Ref. 8). As 
first noted by Isogai (Ref. 5), the flutter mode is essentially the first 
(predominantly bending) natural mode. The flutter frequency is close to the 
first coupled natural frequency Wj/(«) a until the nose of the bend-back is en- 
countered. At this point, the flutter frequency increases to a value between 
the two coupled natural frequencies and the flutter mode also changes, although 
it is still associated with the first predominantly bending branch. 



FIGURE 4 
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TRANSIENT SOLUTIONS FOR CASE A AROUND LOWER FLUTTER BOUNDARY AT N « 0.9 


Typical transient solutions for Case A are shown in Figs. 5,6. At a Mach 
No. of 0.9, multiple flutter solutions occur due to the bend back of the flutter 
boundary (see Fig. 3). Figure 5 illustrates the dynamic behavior of the air- 
foil, plotted as h(t)/b and a(t) vs. time, immediately above and below the 
lower flutter point at M = 0.9 In this case, the airfoil is stable below 
(bottom figure) and unstable above (top figure) the neutral stability boundary 

(CL U„/bw vs. M.) Here, the airfoil has been forced for 3-6 cycles in pure 
F F a 

pitch, with an amplitude of 0.1°, and then released at t=0 . 
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TRANSIENT SOLUTIONS FOR CASE A AROUND UPPER FLUTTER BOUNDARY AT M * 0.9 


In the vicinity of the upper flutter point at M = 0.9, the stability 
behavior is reversed from that observed around the lower point. The airfoil is 
now stable for values of nondimensional airspeed U/b <o a , above the neutral 
stability boundary, as shown in the bottom diagram of Fig. 6. Conversely, the 
airfoil is unstable for values of U/bw a below the flutter boundary. 




TIME 


FIGURE 6 
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FLUTTER CALCULATIONS FOR CASE B 


In Table 2, the results from our present flutter calculations for Case B 
are compared to predictions by previous researchers using various TSD codes. 
This case is the same as Case B considered by Isogai in Ref. 6. Note that the 
present Euler calculations predict a somewhat higher flutter speed than the TSD 
calculation by Isogai, but still below the speed predicted by classical linear 
theory. The flutter speed predicted by Ueda and Dowell (Ref. 9), using the 
describing function method based on LTRAN2 aerodynamics, is significantly below 
the predictions of the Euler code. 


TABLE 2 

Comparison of Predicted Flutter Speed for Case B 


Method 

a 

Up/b(*) a 

2kp 

Present 

0.1° 

3.43 

0.201 

Ueda & Dowell 9 

0.25° (4>i ) 

2.95 

0.221 

Isogai 6 

0.1° 

3.25 

0.215 

Linear Theory 

- 

3 . 86 

0.210 
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NONLINEAR DEPENDENCE OF FLUTTER SPEED 
ON INITIAL FORCING AMPLITUDE FOR CASE C 

In Refs. 9-10, Ueda and Dowell investigated the nonlinear amplitude 
dependence of the flutter boundary for Case C, Table l. They found a distinct 
drop in the flutter speed as the amplitude of the effective induced angle of 

attack <|>i = a + h c /U exceeded about 1°, where h c is the plunging velocity at 
midchord. Figure 7 shows results from the present Euler calculations, plotted 
as flutter speed vs. initial forcing amplitude in pitch (prior to release). 
Note that the flutter boundary is not very sensitive to a in the range 0°-5°, 
and that the results obtained are sensitive to the initial forcing frequency. 




FIGURE 7 
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TRANSIENT SOLUTIONS FOR CASE C 


Typical stable and unstable transient, solutions are shown in Figs. 8 and 9, 
corresponding to initial forcing amplitudes of a = 1° and 4°, respectively. 
The flutter mode is again a predominantly bending mode and emerges quickly 
(within a couple of cycles), despite the fact that the initial disturbance is 
pure torsion. This rapid convergence toward the significant aeroelastic mode 
was also observed in most of the transient solution of Cases A and B as well. 
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FIGURE 8 


FIGURE 9 
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CONCLUSIONS 


1. Typical section flutter calculations based on the two-dimensional unsteady 
Euler equations are now feasible. 

2. Flutter speeds predicted by the present Euler code are in good overall 
agreement with previous TSD calculations, except in cases where strong 
shocks are present. 

3. The Euler code calculations predict a transonic dip similar to the corre- 
sponding dips predicted by TSD codes, but shifted toward higher Mach numbers. 

4. Multiple flutter points occur at certain Mach numbers, caused by a bend- 
back of the flutter boundary. 

5. The amplitude dependence of Up appears to be less than might be expected. 
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